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Abstract. We propose a useful method for mapping large-scale velocity fields in the solar photosphere. It is based on the 
local correlation tracking algorithm when tracing supergranules in full-disc dopplergrams. The method was developed using 
synthetic data. The data processing the data are transformed during the data processing into a suitable coordinate system, the 
noise is removed, and finally the velocity field is calculated. Resulting velocities are compared with the model velocities and 
the calibration is done. From our results it becomes clear that this method could be applied to full-disc dopplergrams acquired 
by the Michelson Doppler Imager (MDI) onboard the Solar and Heliospheric Observatory (SoHO). 
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1. Introduction 

The solar photosphere is a very dynamic layer of the solar 
atmosphere. It is strongly influenced by the underlying con- 
vection zone. Despite years of intensive studies, the velocity 
fields in the solar photosphere remain not very well known. 
The evidence of the vigorous and sometimes chaotic charac- 
ter of the motions of observed structures in the photosphere 
(sunspots, granules, and other features) came already from the 
first systematic studies made in 19th century, of which let us 
at least mention the discovery of the solar diff'erential rotation 
(Carrington I1859> . Motions in the photosphere are strongly 
coupled with magnetic fields. The large-scale velocity fields 
are very important for studies of the global solar dynamo. 

An attempt to describe the differential rotation by 
a parabolic dependence did not make for clear results. 
Coefficients of the parabola differed according to traced objects 
and also changed with time, when tracing one type of object 
(reviewed e. g. by Schroter |l985> . The character of the differ- 
ential rotation has never been in doubt. 

It follows from these arguments that a temporally vari- 
able streaming of the plasma exists on the surface of the Sun, 
which can be roughly described by the differential rotation. 
This streaming has a large-scale character, large-scale plasma 
motions were studied for example on the basis of tracking the 
magnetic structures (e. g. Ambroz |200 1 al 1200 1 b> . The long- 
term Doppler measurements done by the MDI onboard SoHO 
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make it possible to extend the studies of large-scale velocities 
in the solar photosphere. The knowledge of the behaviour of 
velocities in various periods of the solar activity cycle could 
contribute to understanding the coupling between the velocity 
and magnetic fields and of the solar dynamo function. 

There are at least three methods calculating photospheric 
velocities: 



Direct Doppler measurement - provides only one compo- 
nent (line-of-sight) of the velocity vector. These velocities 
are generated by local photospheric structures, amplitudes 
of which are significantly greater than amplitudes of the 
large-scale velocities. The complex topology of such struc- 
tures complicates an utilisation for our purpose. Analysing 
this component in different parts of the solar disc led to 
very important discoveries (e. g. supergranulation - Hart 
[Tgjel and Leighton et al. [T9^ . 

Tracer-type measurement - provides two components of 
the velocity vector. When tracing some photospheric trac- 
ers, we can compute the local horizontal velocity vectors in 
the solar photosphere. Tracking motions of sunspots across 
the solar disc led to the discovery of the differential rotation 
(Carrington[l859). 

Local helioseismology - provides a full velocity vector. 
Although local helioseismology (e. g. Zhao l200'4t is a very 
promising method, it is still in progress and until now does 
not provide enough reliable results. 
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Fig. 1. Comparison of the real dopplergram (left, observed by MDI/SoHO) and the simulated one (right). Both images are visually 
very similar The black colour means line-of-sight velocity -700 m s"' (towards an observer) while white represents +700 m s"' . 



Since the photosphere is a very thin layer (0.04 % of the 
solar radius), the large-scale photospheric velocity fields have 
to be almost horizontal. Then, the tracer-type measurement 
should be sufficient for mapping the behaviour of such veloci- 
ties. In this field the local correlation tracking (LCT) method is 
very useful. 

This method was originally designed for the removal of 
the seeing-induced distortions in image sequences (November 
I1986> and later used for mapping the motions of granules in the 
series of white-light images (November & Simon [T988> . The 
method works on the principle of the best match of two frames 
that record the tracked structures at two different instants. For 
each pixel in the first frame, a small correlation window is cho- 
sen and is compared with a somewhat displaced window of the 
same size in the second frame. A vector of displacement is then 
defined as a difference in the coordinates of the centres of both 
windows when the best match is found. The velocity vector is 
calculated from this displacement and the time lag between two 
frames. 

LCT was recently used for tracking many features in vari- 
ous types of observations, especially for tracking the granules 
in high-resolution white-light images (e. g. Sobotkaetal. [T999l 
I2OOOI . 

The method needs a tracer - a significant structure recorded 
in both frames, the lifetime of which is much longer than the 
time lag between the correlated frames. We decided to use 
the supergranulation pattern in the full-disc dopplergrams, ac- 
quired by the MDI onboard SoHO. We assume that supergran- 
ules are carried as objects by the large-scale velocity field. This 
velocity field is probably located beneath the photosphere, so 
that the resulting velocities will describe the dynamics in both 
the photospheric and subphotospheric layer. The existence of 



the supergranulation on almost the whole solar disc (in contrast 
to magnetic structures) and its large temporal stability make the 
supergranulation an excellent tracer 

The resulting velocities cover the whole solar disc. Hence 
the velocity field also describes a motion of the supergranu- 
lation in the areas occupied by active regions or by magnetic 
field concentrations. This fact allows the results to be used to 
study the mutual motions of substructures like sunspots, mag- 
netic field in active regions, background fields, and the quiet 
photosphere. 

Supergranules are structures with a strong convection cou- 
pling. The mean size of the supergranular cell is approx. 30 Mm 
(e. g. Wang & Zirin I1989t . with the size and shape depen- 
dent on the phase of the solar cycle (e. g. Sykora I1970> . 
Supergranules are quite stable with a mean lifetime of approx. 
20 hours (e. g. Leighton ll964t . Information about the distri- 
bution function of the lifetime was published by DeRosa et 
al. J2000> . The internal velocity field in the supergranular cell 
is predominantly horizontal (e. g. Hathaway et al. I2002> with 
the amplitude approx. 300 ms Due to the horizontality of 
the internal velocity field, the supergranules can be observed in 
dopplergrams on the whole solar disc except for its centre. 

We do not propose that this method capable of measuring 
velocities of order 1 ms"', but we do expect that the large- 
scale velocities will have magnitudes at least one order greater 
We also have to take the largest-scale velocities into account 
like the differential rotation or meridional circulation, which 
have velocities of at least 10 ms"'. If for example we take 
differential rotation described by the formula a;[deg/day] - 
13.064 - 1.69 sin- /7 - 2.35 sin''/? (Snodgrass [T954b . we have 
to expect a velocity approx. 190 ms"' in b = 60° in the 
Carrington's coordinate system. The main goal of our study 
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(and the proposed method of mapping the velocities should be 
proxy for this purpose) is to separate the superposed velocity 
field into the components and to investigate their physics. 

This paper is the first one of a series about the large-scale 
photospheric velocity fields. In the following papers we shall 
apply the method here suggested to observed data and describe 
the properties of the large-scale velocity fields in different pe- 
riods of the solar activity cycle. 

2. Synthetic data 

A recent experience with applying this method to observed 
data (e. g. Svanda et al. 12005 II has shown that for the proper 
setting of the parameters and for the tuning of the method, 
synthetic (model) data with known properties are needed. The 
synthetic data for the analysis come from a simple simulation 
(SISOID code = S'/mulated S'upergranulation as Observed In 
Dopplergrams) with the help of which we can reproduce the 
supergranulation pattern in full-disc dopplergrams. 

The SISOID code is not based on physical principles taking 
place in the origin and evolution of supergranulation, but in- 
stead on a reproduction of known parameters that describe the 
supergranulation. Individual synthetic supergranules are char- 
acterised as centrally symmetric features described by their po- 
sition, lifetime (randomly selected according to the measured 
distribution function of the supergranular lifetime - DeRosa et 
al.|2000^, maximal diameter (randomly according to its distri- 
bution function - Wang & Zirin il989< and characteristic values 
of their internal horizontal and vertical velocity components 
(randomly according to their distribution function - Hathaway 
etal. l200^ . 

The most important simplification in the SISOID code is 
that individual supergranules do not influence each other, but 
simply overlap. The final line-of-sight velocity at a certain 
point is given by the sum of line-of-sight velocities of indi- 
vidual synthetic supergranules at the same position. 

New supergranules can arise inside the triangle of neigh- 
bouring supergranules (identification of such triangles is done 
by the Delaunay triangulation algorithmed by Barrv ll99Tt only 
when the triangle is not fully covered by other supergranules 
and when anyone of the supergranules located at the vertices of 
the triangle is not too young, so that in the future it could fully 
cover the triangle. The position of the origin of the new super- 
granule is the centroid of the triangle; each vertex is weighted 
by the size of its supergranule. 

The SISOID simulation is done in the pseudocylindri- 
cal Sanson-Flamsteed coordinate grid (Calabretta & Greisen 
I2002> : the transformation from the heliographic coordinates is 
given by the formula: 

X — j?cosi^, y — (fi, (1) 

where x and y are coordinates in the Sanson-Flamsteed coor- 
dinate system, and § and tp are heliographic coordinates origi- 
nating at the centre of the disc. At each step an appropriate part 
of the simulated supergranular field is transformed into helio- 
graphic coordinates. The output of the program is a synthetic 
dopplergram of the solar hemisphere in the orthographical pro- 
jection to the disc. We assume in our simulation that the Sun 



lies in an infinite distance from the observer and that the P (po- 
sition angle of solar rotation axis) and bo (heliograhic latitude 
of centre of solar disc) angles are known. 

Each step in calculation includes the evaluation of the pa- 
rameters of individual supergranules, then small old supergran- 
ules under the threshold (2 Mm in size) are removed from the 
simulation and all the triangles are checked, whether a new su- 
pergranule can arise inside them. This step in the SISOID code 
corresponds to 5 minutes in real solar time. The computation is 
always started from the regular grid. Properties of "supergran- 
ules" are chosen randomly according to their real distribution 
functions. The first 1000 steps are "dummy", i. e., no vector 
velocity field is included and no synthetic dopplergram is cal- 
culated. This starting interval is taken for the stabilisation of the 
supergranular pattern. In the next steps, the model vector veloc- 
ity field has already been introduced. This field influences only 
the positions of individual cells. The dopplergram is calculated 
every third step. For one day in real solar time, 96 dopplergrams 
are calculated. 

The model velocity field with Carrington rotation added is 
applied according to the assumption of the velocity analysis 
that the supergranules are carried by a velocity field on a larger 
scale. In the simulation only, the position of individual super- 
granules is influenced, and no other phenomena are taken into 
account. These synthetic dopplergrams are visually similar to 
the real observed dopplergrams (see Fig.Q. 

3. Method of data processing 

The MDI onboard SoHO acquired the full-disc dopplergrams at 
a high cadence in certain periods of its operation - one obser- 
vation per minute. These campaigns were originally designed 
for studying the high-frequency oscillations. The primary data 
contain lots of disturbing effects that have to be removed be- 
fore ongoing processing; the rotation line-of-sight profile, p- 
modes of solar oscillations. We detected some instrumental ef- 
fects connected to the data-tranfer errors. It is also known that 
the calibration of the MDI dopplergrams is not optimal (e. g. 
Hathaway et al. l2002> and has to be corrected to avoid system- 
atic errors. While examining long-term series of MDI dopp- 
lergrams, we have met systematic errors connected to the re- 
tuning of the interferometer. We should also take those geo- 
metrical effects into account (finite observing distance of the 
Sun, etc.) causing bias in velocity determination. According 
to Strous J2000> for example, the bias coming out of a per- 
spective is about 2 ms and it depends on the position on 
the disc. It has been proven (e. g. Liu & Norton I200H that 
MDI provides reliable velocity measurements when the mag- 
netic field is lower than 2000 Gauss. The velocity observation 
by MDI will induce up to 100 % error if the magnetic field is 
higher than 3000 Gauss due to the magnetic sensitivity of the 
used Ni I line and the limitations of computational algorithm, 
which cause crosstalks between measured MDI dopplergrams 
and magnetograms. The removal of these effects will be de- 
scribed in detail in the next paper, while the synthetic data used 
in this study do not suffer from these phenomena. 

As input to the data processing we take a one-day observa- 
tion that contains 96 full-disc dopplergrams in 15-minute sam- 
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Fig. 2. Left - model vector velocity field. Right - a velocity field that was computed by applying the LCT method to the synthetic 
supergranulation pattern with the imposed model field. The arrow lengths, representing the velocity magnitudes, have the same 
scale. The images are visually very similar, however the magnitudes of calculated velocities are underestimated. 



pling. Structures in these dopplergrams are shifted with respect 
to each other by the rotation of the Sun and by the velocity field 
under study. 

First, the shift caused by the rotation has to be removed. For 
this reason, the whole data series (96 frames) is "derotated" 
using Carrington's rotation rate, so that the heliographic lon- 
gitude of the central meridian is equal in all frames and also 
equals the heliographic longitude of the central meridian of the 
central frame of the series. This data-processing step causes the 
central disc area ("blind spot" caused by prevailing horizontal 
velocity component in supergranules) in the derotated series to 
move with the Carrington's rate. During the "derotation" the 
seasonal tilt of the rotation axis towards the observer (given 
by b() - heliographic latitude of the centre of the disc) is also 
removed, so that bo = in all frames. 

Then the data series is transformed to the Sanson- 
Flamsteed coordinate system to remove the geometrical dis- 
tortions caused by the projection of the sphere to the disc. 
Parallels in the Sanson-Flamsteedpseudocylindrical coordinate 
system are equispaced and projected at their true length, which 
makes it an equal area projection. Formulae of the transforma- 
tion from heliographic coordinates are given by Eq. Q. 

The noise coming from the evolutionary changes in the 
shape of individual supergranules and the motion of the "blind 
spot" in the data series with the Carrington's rotational rate 
are suppressed by the k-cj filter in the Fourier domain (Title 
et al. 119891 Hirzberger et al. I1997> . The cut-off velocity is set 
to 1 500 ms ' and has been chosen on the basis of empirical 
experience. 

The existence of the differential rotation complicates the 
tracking of the large-scale velocity field, because the ampli- 



tudes and directions of velocities of the processed velocity field 
have a significant dispersion. We have found that, when the 
scatter of magnitudes is too large, velocities of several hun- 
dred m s ' cannot be measured precisely by the LCT algorithm 
where the displacement limit for correlation was set to detect 
velocities of several tens of m s ' . Therefore the final velocities 
are computed in two steps. The first step provides a rough in- 
formation about the average zonal flows using the differential 
rotation curve 

Lo^A+B sin- b + C sin'* b , (2) 

and calculating its coefficients. 

In the second step this average zonal flow is removed from 
the data series, so that during the "derotation" of the whole 
series the differential rotation inferred in the first step and ex- 
pressed by (13 is used instead of the Carrington rotation. The 
scatter of the magnitudes of the motions of supergranules in the 
data transformed this way is much smaller, and a more sensitive 
and precise tracking procedure can be used. 

The LCT method is used in both steps. In the first step, the 
checked range of velocity magnitudes is set to 200 ms"', but 
the accuracy of the calculated velocities is roughly 40 ms"'. 
In the second step the range is only 100 ms"' with much bet- 
ter accuracy. The lag between coiTelated frames equals in both 
cases 16 frame intervals (i. e. 4 hours in solar time), and the 
correlation window with FWHM 30 pixels equals 60" on the 
solar disc in the linear scale. In one observational day, 80 pairs 
of velocity maps are calculated and averaged. 

For the calculation we use the adapted program 
flowmaker .pro originally written in IDL by Molowny- 
Horas & Yi (I1994> . The algorithm has a limitation in the range 
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of displacements that are checked for each pixel. The quality of 
correspondence (in our case the sum of absolute differences of 
both correlation windows) is computed in nine discrete points, 
then the biquadratic surface is fitted through these nine points, 
and an extremum position (Darvann I1991I I is calculated. The 
final displacement vector is equivalent to the position of the 
extremum. 

4. Results 

In our tests we have used lots of variations of simple axisym- 
metric model flows (with a wide range of values of parame- 
ters describing the differential rotation and meridional circula- 
tion) with good success in reproducing the models. When com- 
paring the resulting vectors of motions with the model ones, 
we found a systematic offset in the zonal component equal to 
Voffset, zonal = - 15 m s"' . This constant offset appeared in all the 
tested model velocity fields and comes from the numerical er- 
rors during the "derotation" of the whole time series. For the 
final testing, we used one of the velocity fields obtained in our 
previous work (Svanda et al. l2005t . This field approximates the 
velocity distribution that we may expect to observe on the Sun. 
The model flows have structures with a typical size of 60", 
since they were obtained with the correlation window of this 
size. 

The calculated velocities (with Vofjset, zonal - -15 ms"' 
corrected) were compared with the model velocities (Fig. |2}. 
Already from the visual impression it becomes clear that 
most of vectors are reproduced very well in the direction, 
but the magnitudes of the vectors are not reproduced so well. 
Moreover, it seems that the magnitudes of vectors are underes- 
timated. This observation is confirmed when plotting the mag- 
nitudes of the model vectors versus the magnitudes of the cal- 
culated vectors (Fig. |3}- The scatter plot contains more than 
1 million points, and most of the points concentrate along a 
strong linear dependence, which is clearly visible. This depen- 
dence can be fitted by a straight line that can be used to derive 
the calibration curve of the magnitude of calculated velocity 
vectors. The calibration curve is given by the formula 

Vcoi- = 1.13 Vcalc, (3) 

where Vcaic is the magnitude of velocities coming from the LCT, 
and Vcoi the corrected magnitude. The directions of the vectors 
before and after the correction are the same. The uncertainty of 
the fit can be described by 1-cr-error 15 ms"' for the velocity 
magnitudes under 100 ms ' and 25 ms"' for velocity mag- 
nitudes greater than 100 ms '. The uncertainties of approx. 
15 ms"' have their main origin in the evolution of supergran- 
ules. We studied the dependence of the error of velocity deter- 
mination on the time lag used when no model velocity field was 
introduced. We found that this dependence is slowly increasing 
with the time lag (Fig. ^ due to the evolution of individual 
supergranules. Evolution of supergranules is only one part of 
story, but it gives the lower limit of accuracy that can be ob- 
tained by this method. We also ran a test of the LCT sensitivity 
on the evolution of supergranules when a known underlying 
velocity field is introduced and came to similar results. 
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Fig. 3. Scatter plot for inference of the calibration curve. 
Magnitudes of calculated velocities are slightly underestimated 
by LCT, but the linear behaviour is clearly visible. A line rep- 
resenting the 1 : 1 ratio is displayed. The calibration affects only 
the magnitudes of the flows, while the directions do not need 
any correction. 
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Fig. 4. Dependence of the 1-cr-error of the calculated velocity 
on the time lag when no model velocity field was introduced 
and only the evolution of supergranules have been taken into 
account. The dashed line represents lag 16 (4 hours) that is usu- 
ally used in our method. 

We tested the sensitivity of the method to the choice of val- 
ues of FWHM of the correlation window and of the lag be- 
tween correlated frames in the LCT method. We found that our 
method is practically insensitive to the choice of the time lag 
between correlated frames when the lag is in the interval of 10- 
24 (2.5-6 hours). The larger the lag we choose, the lower the 
velocities we are able to detect. On the other hand, we have to 
take into account that a larger lag between correlated frames 
causes more noise in calculated results coming from evolution- 
ary changes of supergranules and probably also from evolution- 
ary changes in the velocity field under study. According to our 
tests, for a time lag greater than 30 (7.5 hours), the numerical 
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Fig. 5. Influence of the calculated velocity field on the choice of FWHM of the correlation window. Left 
The model velocity field is the same as in Fig.|2 



120", right - 200". 



noise raises very fast. The lag 16 (4 hours) seems to be a good 
tradeoff^ between sensitivity and noise. 

The choice of different FWHMs of the correlation window 
changes the spatial resolution according to FWHM. The gen- 
eral character of the vector field is preserved within the limits of 
resolution (cf. Fig.|5}- Larger FWHM causes a smoothing of re- 
sults and an underestimation of vector magnitudes (cf. Fig.|6}. 
We found the used parameters (FWHM 60", lag 4 hours) to 
be the best compromise, however these values can be changed 
during the work on real data. 
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Fig. 6. 

of the 



Histograms of velocity magnitudes for various FWHM 
LCT algorithm. 



5. Conclusions 

We have developed a method for mapping the velocity fields 
in the solar photosphere based on the local correlation track- 
ing (LCT) algorithm. The method consists in two main steps. 
In the first main step the mean zonal velocities are calculated 
and, on the basis of expansion to the Fay's formula (|2}, the 
differential rotation is removed. In the second main step, the 
LCT algorithm with an enhanced sensitivity is applied. Finally, 
the differential rotation (obtained in the first step) is added to 
the vector velocity field obtained in the second main step. Both 
main steps can be divided into a few substeps, which are mostly 
common. 

1. "Derotation" of the data series using the Carrington's rota- 
tion rate in the first step and using the calculated differential 
rotation in the second step. 

2. Transformation of the data series into the Sanson- 
Flamsteed coordinate system to remove the geometrical 
distortion caused by the projection to the disc. 

3. The k-u) filtering (with the cut-off velocity 1500 ms ') 
for suppression of the noise coming from the evolutionary 
changes of supergranules, of the numerical noise, and for 
the partial removal of the "blind spot" (an effect at the cen- 
tre of the disc). 

4. Application of the LCT; the lag between correlated frames 
is 4 hours, the correlation window with FWHM 60", the 
measure of correlation is the sum of absolute differences 
and the nine-point method for calculation of the subpixel 
value of displacement is used. The calculated velocity field 
is averaged over the period of one day. 

Finally, the magnitude of the calculated vector field is corrected 
using formula (jSjl obtained from the tests on the synthetic data. 
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According to our test, the method provides a very reliable 
tool for mapping the velocity fields with the spatial resolution 
of 60" (43,5 Mm) on the solar disc and with the accuracy of 
15 ms"' for velocity magnitudes under 100 m s"' and 25 m s"' 
for velocity magnitudes greater than 100 m s ' . The method is 
ready to use on the real full-disc dopplergrams acquired by the 
MDI onboard SoHO. Our method shows that the usability of 
the LCT method is wider than originally expected. 
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